Normal Stress Distribution of Rough Surfaces in 
Contact 

Alex Hansen, 1 ' 2 Jean Schmittbuhl, 3 G. George Batrouni, 4 and Fernando A. 
de Oliveira 

International Center for Condensed Matter Physics, Universidade de Brasilia, Brazil 

Abstract. We study numerically the stress distribu- 
tion on the interface between two thick elastic media 
bounded by interfaces that include spatially correlated 
asperities. The interface roughness is described using 
the self-affine topography that is observed over a very 
wide range of scales from fractures to faults. We analyse 
the correlation properties of the normal stress distribu- 
tion when the rough surfaces have been brought into 
full contact. The self affinity of the rough surfaces is 
described by a Hurst exponent H . We find that the 
normal stress field is also self affine, but with a Hurst 
exponent H — 1. Fluctations of the normal stress are 
shown to be important, especially at local scales with 
anti-persistent correlations. 



Theories describing the elastic properties of two me- 
dia in contact through rough surfaces have important 
applications in a wide range of geophysical problems, 
such as earthquakes, fracture, fluid permeability or rock 
friction [Scholz, 1990]. Asperities exist at all scales: 
grain roughness is relevant for closure of rock joints 
[Brown and Scholz, 1986] and seamounts might induce 
large scale stress fluctuations along subduction slabs 
[Dmowska et at, 1996]. Whatever the scale of the asper- 
ities in contact is, when they are attached to an elastic 
medium and are loaded, they interact and concentrate 
high stresses. Friction properties of an interface are 
very dependent on the heterogeneities of the normal 
stresses [Dieterich and Kilgore, 1996]. At fault scale, 
residual stresses resulting from asperity squeeze might 
be responsible for heterogeneities of the dynamic stress 
field and influence the earthquake propagation [Bou- 
chon et at, 1998]. 

We consider in this letter the normal component of 
the stress field ctn that appears on the interface between 
two elastic media with rough surfaces (see Figure 1) 
[Sayles, 1996]. We assume that possible local plastic 
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deformations where the elastic yield stress of asperities 
is overcome, have a negligible effect on the stress distri- 
bution along the interface. The roughness is assumed 
to be self affine. The surface is described as h(x,y). 
Rescaling the two coordinates x — > Ax and y — ► Ay, 
necessitates a rescaling of the height h — ► A H [Feder, 
1988]. The surface is then self affine with a Hurst ex- 
ponent H . In Figure 1 an example of surfaces with this 
property is shown. 

This surface was obtained from an impact fractur- 
ing of a granite block (25cm x 25cm) [Lopez and Schmit- 
tbuhl, 1998]. A wide range of experimental data [Brown 
and Scholz, 1985a; Power et al, 1987; Bouchaud et al, 
1990; Mdl0y et al, 1992; Schmittbuhl et al, 1995] sup- 
port the hypothesis that not only are surfaces produced 
by brittle fracture self affine, but their Hurst exponent 
generally equals H = 0.80 independently of the mate- 
rial [Bouchaud, 1997]. We have measured the roughness 
of the granite block using a profilometer and analysed 
the spatial correlations of the surface with the aver- 
age wavelet coefficient technique [Mehrabi et al, 1997; 
Simonsen et al, 1998]. This consists in wavelet trans- 
forming each one-dimensional trace h(x, y = const) us- 
ing the Daubechie-12 wavelet basis and averaging over 
the wavelet coefficients at each length scale b — 2 k , 
where k is an integer. If the trace is self affine, the 
averaged wavelet coefficients scale as 



Ah ~ b H+1 ' 2 



(1) 



where Ah is an average over the position of the wavelet. 

We show in Figure 2 the average wavelet coefficients 
vs. b for the granite surface of Figure 1. The slope of the 
least-squares fit is 1.30 ± 0.02, giving H = 0.80 ± 0.02. 

Integrating the Lame equations for an infinite block 
limited by an infinite plane (x, y, z = 0), gives the Green 
function G for the deformational response u in the z 
(vertical) direction at a point (x u , y u ) in the plane z = 
from a distribution f(xf,yf) of applied forces in the 
vertical direction: 



u(x u ,y u )= J J G(x u -x f ,y u -y f )f(x f ,y f )dx f dy f 
where [Landau and Lifshitz, 1958] 



G(x u -x f ,y u - y f ) 



l-s 2 1 



nE r 



(2) 
(3) 



Here r = (x u — Xf) 2 + (y u — yf) 2 , E is the elastic 
constant and s is the Poisson ratio. Deformation also 
occurs within the (x, y)-plane when vertical force is ap- 
plied. However, these fall of as 1/r 2 . Consequently, we 
ignore them compared to the deformation in the vertical 
direction. 

We are interested in rough self-affine surfaces. How- 
ever, with a Hurst exponent H < 1, the surfaces are 
asymptotically flat. This can be easily seen by calcu- 
lating the rms fluctuations of the surface, w 2 = {{h — 
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{h)) 2 )LxL, where (• • -)lxl is an average over an area of 
size L x L of the plane (x,y). When the surface is self 
afhne, we have w <~ L H leading to w/L ~ \/L}~ H — > 
as L — > oo. Thus, it is asymptotically flat. As we are 
primarily interested in the scaling properties of the nor- 
mal component of the stress field o~n on large scales, it 
is a good approximation to use the flat-surface Green 
function, (3). This approximation also leads to the force 
component in the vertical direction being directly pro- 
portional to the normal stress. 

The problem we have set out to study is that of two 
self-affine rough surfaces in full contact. However, as- 
suming that one of the surfaces is rough and infinitely 
hard, and the other clastic and initially flat, we find 
the same normal stress field at the interface as in the 
original problem within the approximation using the 
flat-surface Green function, (3) and using the compos- 
ite topography introduced by Brown and Scholz [1985b] 
(i.e. the sum of both topographies). We will, therefore, 
study this second problem since it is easier to implement 
numerically 

When the two surfaces are in full contact, the defor- 
mation field u will be equal to minus the local height, 
u = —h, when in-plane deformations are ignored. Thus, 
the deformation field is self-affine, with a Hurst expo- 
nent H . 

Since Eq. (2) is linear, simple scaling arguments tell 
us how the force field / scales given that u is self affine. 1 
If we scale {x u ,y u ) — > (Xx u , Xy u ) and (x fl y f ) — > (Xxf,Xy f ), 
Eqs. (2) and (3) immediately gives the scaling relations 



W) 4ttE 



Ax -> XAx , 
Ay -» XAy , 
u — > X H u , 
G — > A _1 G , 



(4) 



with Ax = x u — Xf and Ay = y u — Uf- Thus, the 
force field / is self affine with a Hurst exponent equal 
to H a = H-l. 

In order to demonstrate the validity of Eq. (4), we 
solve Eq. (2) numerically for /. This is done in Fourier 
space (using FFTs [Stanley and Kato, 1997]) since the 
Green function is diagonal there. We start out by 
defining the Green function on the L x L square lat- 
tice as follows. For each node (i,j), we define r± — 
((i - l) 2 + (j - l) 2 ) 1 / 2 , r 2 = ((i - l) 2 + (L + 1 - 3 ?) XI \ 
r 3 = ((L + l-if + ij-l) 2 ) 1 ' 2 , andr 4 - {{L + l-i) 2 + 
(L + 1 — j) 2 ) 1 / 2 . The Green function on this lattice is 
then given by 



1 111 

1 7 + - + - + - 

max(ri , e) ri rz 



■ (5) 



1 lf the underlying equations were not linear, much more power- 
ful methods to determine the scaling behavior would be necessary 
such as functional renormalization [Barabdsi and Stanley, 1997]. 



Thus, the singularity of the Green function is situated 
at [i = l,j = 1). We have introduced a cutoff in r\ 
equal to e. We choose it to be a quarter of the lattice 
spacing, i.e., 1/4. The reason for introducing the three 
other radii r 2 , r 3 , and is that the Fourier transform 
makes the lattice periodic. The three additional radii 
signify the mirror image of the singularity resulting from 
one reflection — we do not introduce further reflections 
since, with our choice of parameters, their effect is neg- 
ligible. The deformation field u and the Green fuction, 
G were then Fourier transformed, and Eq. (2) solved in 
Fourier space. The resulting force field was then Fourier 
transformed back to real space. To within the approxi- 
mations we have made, the force field is proportional to 
the normal stress field <7jv ■ We show in Figure 3 the nor- 
mal stress field corresponding to the full contact of the 
fracture in Figure 1. The strong small-scale variations 
in the normal stress distribution are consistent with the 
observations reported in Mendelsohn et al. [1998]. 

In Figure 2, we show the wavelet analysis of the nor- 
mal stress field obtained for the granite surface. The 
least-squares fit gives a slope of 0.30 ± 0.04 correspond- 
ing to a Hurst exponent of —0.20 ± 0.04. Thus, the 
relation H a = H — 1 is supported. 

In order to study systematically the relation between 
the Hurst exponent of the deformation field and that of 
the normal stress field, we have generated artificial self- 
affine surfaces, using the Fourier method [Sahimi, 1998]. 
This allows us to generate and subsequently average our 
results over many surfaces for each Hurst exponent, in 
practice 100 surfaces. After obtaining the stress fields 
for each surface, we analysed spatial correlations of both 
the surfaces and the stress fields with two techniques: 
the average wavelet spectrum of one-dimensional traces 
obtained by cutting the surface along lines and the two 
dimensional Fourier spectrum of the full surface. Fig- 
ure 2 presents the average wavelet spectra of the sur- 
faces and their corresponding stress field for the syn- 
thetic fracture surfaces with the same Hurst exponent as 
the one that is observed for fracture surfaces: H = 0.80. 
We treated each surface and its corresponding stress 
field as consisting of 2048 one-dimensional traces, and as 
there were 100 surfaces, our averages are over 100 x 2048 
one-dimenensional traces. The scaling of the synthetic 
surfaces is in good agreement with that of the mea- 
sured surface. Computed full contact stress fields of 
both types of surfaces are also in good agreement sup- 
porting the relation: H a = H — 1. The two dimensional 
Fourier spectrum is computed from the two dimensional 
Fourier transform of the surface and is expected to scale 
for self-affine surfaces as [Sahimi, 1998]: 

H.I-2-2H 



P(\k\ 



(6) 



We generalize the analysis for different Hurst ex- 
ponents that describe different spatial correlations be- 
tween asperities. In Figure 4, we show H a as a function 
of H for the artificially generated 2048 x 2048 surfaces 
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and analyses with both techniques. The straight line 
corresponds to 

H„ = H-1. (7) 

We see that the numerical results and Eq. (7) are in 
excellent agreement for the two dimensional Fourier 
spectrum. The agreement is good with the one dimen- 
sional technique only for sufficiently large roughness ex- 
ponents of the rough surface. We emphasize that mea- 
surements of low roughness exponent have to be done 
with two dimensional techniques [Hansen et al., 2000]. 

The Hurst exponent is directly related to the spa- 
tial correlations of the surface. Eq. 6 shows that for: 
H = — 1 surfaces have a flat spectrum that is are 
white noise with no spatial correlations of the asperities. 
When the Hurst exponent H increases, relative magni- 
tudes of low frequency modes also increase. Asperities 
are smoother and the surface roughness appears more 
and more correlated at large scales. Equation (7) shows 
that the stress field can be calculated approximately as 
a simple derivative of the deformation field. Also, fluc- 
tations of the stress field are significantly higher than 
the deformation field. 
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Figure Captions 



Figure 1. A rough granite surface produced by cleav- 
age. The size is 512 x 512 points. Physical dimensions 
are 12.8mm x 12.8mm. 



Figure 2. The average wavelet coefficient At, for the 
fracture surface of Figure 1 as a function of length scale 
b is shown as filled circles. Empty circles correpond to 
the average over 100 synthetic surfaces of 2048x2048. 
The slope of the straight solid line is 1.30, correspond- 
ing to a Hurst exponent H = 0.80. Average wavelet 
coefficient At, for the stress field shown in Figure 3 as a 
function of length scale b using filled diamonds. Empty 
diamonds describe the analysis of the squeeze of the syn- 
thetic surfaces. The slope of the dashed straight line is 
0.30, corresponding to a Hurst exponent H a = —0.2. 

Figure 3. The numerically calculated stress field on 
the interface between the granite block of Figure 1 when 
brought in complete contact with an infinitely hard flat 
surface. 

Figure 4. Hurst exponent of stress field, H a as func- 
tion of Hurst exponent of deformation field H . The data 
are based on averaging over several samples of artifially 
generated rough surfaces (o for ID wavelet spectrum 
and o for 2D Fourier spectrum). The straight line is 
H a =H-l. 
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Figure 1. A rough granite surface produced by cleavage. The size is 512 x 512 points. Physical 
dimensions are 12.8mm x 12.8mm. 
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Figure 2. The average wavelet coefficient for the fracture surface of Figure 1 as a function of 
length scale b is shown as filled circles. Empty circles correpond to the average over 100 synthetic 
surfaces of 2048x2048. The slope of the straight solid line is 1.30, corresponding to a Hurst 
exponent H = 0.80. Average wavelet coefficient for the stress field shown in Figure 3 as a 
function of length scale b using filled diamonds. Empty diamonds describe the analysis of the 
squeeze of the synthetic surfaces. The slope of the dashed straight line is 0.30, corresponding to 
a Hurst exponent H a = —0.2. 




Figure 4. Hurst exponent of stress field, H a as function of Hurst exponent of deformation field 
H. The data are based on averaging over several samples of artifially generated rough surfaces 
(o for ID wavelet spectrum and o for 2D Fourier spectrum). The straight line is H a = H — 1. 



